
    ##############################################################################################
    # R code to replicate
    # Figure 3 of Measuring Foreign Policy Positions of Members of the US Congress
    ##############################################################################################

    #setwd("C:\\0Jeong\\1_Projects\\Foreign_Policy\\0_Estimating_FP_Ideologies\\2015_PSRM\\Replication\\Replication_Material\\Supplemental_Material")

    Data <- read.table("Figure3_Data1.txt", header=T, sep="") # Roll Call Data
    Time <- read.table("Figure3_Data2.txt", header=T, sep="") # Time Index to be fed into MCMCdynmicIRTi1
    info2 <- read.csv("Figure3_Data3.csv", header=T)          # Legislator Information: ICPSR_ID, congress, party ID, and Name (Adapted from Voteview.com data)
    NOM_Dem.medians <- read.csv("Figure3_Data4.csv", header=T)$Dem.median  # Democratic median score of DW-NOMINATE dimension 1 (computed from Voteview.com data)
    NOM_Rep.medians <- read.csv("Figure3_Data4.csv", header=T)$Rep.median  # Republican median score of DW-NOMINATE dimension 1 (computed from Voteview.com data)


    library(MCMCpack)

    out <- MCMCdynamicIRT1d(Data, 
                            item.time.map=Time$time,
                            mcmc=10000, burnin=1000, thin=10,       
                            verbose=1000, store.item=TRUE,         # I set "store.item=TRUE" to produce outputs for Figure 6 in Supplemental Files. For a quick result, set it FALSE.          
                            theta.constraints=list(A9369A="+",     # Strom Thurmond
                                                   A14105A="+",    # Jesse Helms
                                                   A14920A="-",    # John Kerry 
                                                   A10808A="-"))   # Edward Kennedy

 
    out_ideal <- out[seq(from=1, to=nrow(out), length=min(1000, nrow(out))) , 1:3350]   # select ideal points only
    N <- ncol(out_ideal)

    # to check if converged, select 5 chains randomly.
    postscript("Supplemental_Figure15.eps", height=11, width=7.5, horizontal=F)
    n.size <- 5
    par(mfrow=c(n.size,2))
    for(i in sample(1:3350, n.size, replace=F)){
    plot(density(out_ideal[,i]), main="", xlab="")
    title(main=paste("Ideal Point of", colnames(out_ideal)[i]), font.main=1)
    plot(out_ideal[,i], type="l", xlab="iteration", ylab="")
    title(main="Traceplot", font.main=1)
    }
    dev.off()

    ### Geweke Diagnostics 
    # Note Geweke stat does not work well when the autocorrelation is high. 
    # If so, increase mcmc iterations and the thinning interval (say, from 10 to 50 or 100). (but it means longer simulations).
    # Also, make sure to change "store.item=TRUE" to FALSE to speed up simulation. 

    Geweke.stat1 <- rep(NA, N)
    for(i in 1:N){
    Geweke.stat1[i] <- geweke.diag(out_ideal[,i], frac1=0.4, frac2=0.6)$z
    }
    postscript("Supplemental_Figure16.eps", height=7, width=11, horizontal=F)
    hist(Geweke.stat1, breaks=100, xlim=c(-5, 5), xlab="Geweke Statistics", main="")
    abline(v=c(-1.96, 1.96), lty=2)
    dev.off()

    ##############################################
    ### To save bill parameters for Figure 6 of Supplemental Files
    
    # Store only Difficulty parameters
    out_alpha <- out[seq(from=1, to=nrow(out), length=min(1000, nrow(out))), 3351:7232] 
    dim(out_alpha)
    
    # Store only Discrimination parameters
    out_beta <- out[seq(from=1, to=nrow(out), length=min(1000, nrow(out))), 7233:11114] 
    dim(out_beta)

    mean_alpha <- apply(out_alpha, 2, FUN=mean)
    mean_beta  <- apply(out_beta, 2, FUN=mean)
    
    mean_bill <- cbind(mean_alpha, mean_beta)

    ### Decode the lable into Congress and vote number  
    bill.label <- rownames(mean_bill)    
    congress <- vote.no <- time <- b22 <- b23 <- rep(NA, length(bill.label))
    for(j in 1:length(bill.label)){
    b22[j] <- strsplit(bill.label, split="_v")[[j]][1]
    vote.no[j] <- as.numeric(strsplit(bill.label, split="_v")[[j]][2])
    congress[j] <- as.numeric(strsplit(b22, split="alpha.h")[[j]][2])
    }
    time <- congress - 78
    Vote.Info <- data.frame(bill.label, time, congress, vote.no)
    colnames(Vote.Info) <- c("bill.label","time","congress","vote.no") 
    
    ### Merge Vote Info and the Summary of Bill Parameters
    Bill.Param <- data.frame(Vote.Info, mean_bill)
    write.csv(Bill.Param, "Bill.Param.csv", row.names=F)
    ###################################################  

    ### Prepare the MCMC chains to merge with legislator information
    original.label <- colnames(out_ideal)
    MC_Data <- data.frame(t(out_ideal))
    MC_Data$original.label <- original.label
   
    ### Decode the lable into ICPSR_ID and Congress
    congress <- icpsr_id <- time <- a22 <- a23 <- rep(NA, length(original.label))
    for(j in 1:length(original.label)){
    a22[j] <- strsplit(original.label, split="A.t")[[j]][1]
    time[j] <- as.numeric(strsplit(original.label, split="A.t")[[j]][2])
    icpsr_id[j] <- as.numeric(strsplit(a22, split="theta.A")[[j]][2])
    }
    congress <- time + 78
    
    Info_dynamic <- data.frame(original.label, icpsr_id, time, congress)
    colnames(Info_dynamic) <- c("original.label","ICPSR_ID","time","congress") 
    
    ### Merge with Legislator Information
    Info <- merge(Info_dynamic, info2, by=c("ICPSR_ID", "congress"), all=F)


    ### Save summary outputs to be used for Figures 2 and 4
    summary_results1 <- cbind(apply(out_ideal, 2, FUN=mean), t(apply(out_ideal, 2, FUN=quantile, probs=c(0.025, 0.975))))
    summary_results <- cbind(row.names(summary_results1), summary_results1)
    colnames(summary_results) <- c("original.label", "IRT_mean", "IRT_2.5Q", "IRT_97.5Q")
   
    # Merge Legislator Info and the Summary of Estimates
    Combined <- merge(summary_results, Info, by="original.label", all=F)
    Combined1 <- Combined[order(Combined$congress),]
    write.csv(Combined1, "Summary_Outputs.csv")

    
    ### Merge Legislator Info and the Summary of Estimates
     
    Outputs <- merge(Info, MC_Data, by="original.label", all=F)
    info_columns <- 6 # the number of columns that are NOT MCMC chains

    ### Compute the median of each party
    
    T <- 33 # Number of Congresses
    iter <- nrow(out_ideal) # number of MCMC iterations

    ### Storages
    Dem.median <- matrix(NA, T, iter)
    Rep.median <- matrix(NA, T, iter)

    ### Loop over congresses

    for (t in 1:T){
    Est <- subset(Outputs, time==t)
    
    for (j in 1:iter){
    Dem.median[t,j] <- median(Est[,info_columns + j][Est$party==100])
    Rep.median[t,j] <- median(Est[,info_columns + j][Est$party==200])
    }
    }
    
    ### End of the loop
     
    Dem.median.quant <- apply(Dem.median, 1, FUN=quantile, na.rm=T, probs=c(0.025, 0.5, 0.975))
    Dem.median.mean <- apply(Dem.median, 1, FUN=mean, na.rm=T)
    Dem.median.CI <- rbind(Dem.median.mean, Dem.median.quant)
    
    Rep.median.quant <- apply(Rep.median, 1, FUN=quantile, na.rm=T, probs=c(0.025, 0.5, 0.975))
    Rep.median.mean <- apply(Rep.median, 1, FUN=mean, na.rm=T)
    Rep.median.CI <- rbind(Rep.median.mean, Rep.median.quant)
 

    ### Party medians Comparison with DW-NOMINATE (Computed from the data from Voteview.com)
    year.odd <- seq(from=1945, to=2010, by=2)
        
    postscript("Figure3.eps", height=10, width=8, horizontal=F)
    par(mfrow=c(2, 1))
    plot(1:T, Dem.median.CI[1,], type='n', col='blue', ylim=c(-2,2), xaxt="n", xlab="", ylab="")
    points(1:T, Dem.median.CI[1,],  type='o', lty=1, col='blue', cex=.7)
    segments(1:T, Dem.median.CI[2,], 1:T, Dem.median.CI[4,], col='blue')
    points(1:T, Rep.median.CI[1,],  type='o', lty=2, col='red', cex=.7)
    segments(1:T, Rep.median.CI[2,], 1:T, Rep.median.CI[4,], col='red')
    axis(side=1, las=3, at=1:T, labels=year.odd, cex.axis=.8)
    title(main="Bayesian IRT Estimate",  font.main=1)
    legend('topleft', legend=c("Democrats", "Republicans"), lty=c(1, 2), col=c('blue', 'red'),cex=.8)

    plot(1:T, NOM_Dem.medians, type='o', col='blue', ylim=c(-1,1), xaxt="n", xlab="",ylab="")
    points(1:T, NOM_Rep.medians,  type='o', lty=2, col='red')
    axis(side=1, las=3, at=1:T, labels=year.odd, cex.axis=.8)
    title(main="DW-NOMINATE",  font.main=1)
    legend('topleft', legend=c("Democrats", "Republicans"), lty=c(1, 2), col=c('blue', 'red'),cex=.8)
    dev.off()
       

